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Abstract 

The development of new techniques to quantitatively measure gene 
expression in cells has shed light on a number of systems that display 
oscillations in protein concentration. Here we review the different mech- 
anisms which can produce oscillations in gene expression or protein con- 
centration, using a framework of simple mathematical models. We focus 
on three eukaryotic genetic regulatory networks which show "ultradian" 
oscillations, with time period of the order of hours, and involve, respec- 
tively, proteins important for development (Hesl), apoptosis (p53) and 
immune response (NF-kB). We argue that underlying all three is a com- 
mon design consisting of a negative feedback loop with time delay which 
is responsible for the oscillatory behaviour. 



1 Introduction 

Biological systems display fascinating spatial and temporal patterns, which hint 
that the underlying cellular processes are highly dynamic and operate on a 
wide range of time and length-scales. This is indeed confirmed by measure- 
ments of the temporal dynamics of protein concentrations and gene expression 
levels for various signalling and response systems, made using pulse-labelling 
[T] , /3-galactosidase measurements and immunoblotting [2] , electromobility shift 
assays [3J, real time PCR 0], fluorescence techniques [5J |BJ [71 |S], chromatin 
immunoprecipitation assays [9] and microarrays [10] . Overall, it is now evident 
that regulatory and signal transduction networks do not depend merely on shift- 
ing the relevant protein concentrations from one steady state level to another. 
Rather, the signals often have a significant temporal variation that carries much 
more information and propogates through the regulatory networks in a complex 
manner. With time-resolved data now available for a number of response and 
signalling systems, it is perhaps the appropriate time to explore whether there 
are any commonalities, or "design principles" , in the underlying mechanisms. In 
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this review we will show that a prominent subclass of such systems does indeed 
have a common underlying design structure which combines a negative feedback 
loop with a time delay. 

This subclass consists of systems that display oscillatory behaviour under 
some conditions. The most obvious examples are circadian rhythms and cell 
division. Oscillations are also seen in the levels of cellular calcium [TT]. and 
in embryo development. Somite segmentation, for instance, exhibits clearly pe- 
riodic spatial patterns which are produced by periodic temporal variation of 
proteins like Hesl, Axin, Notch and Wnt [TH [OH Q31 [T5] . We also include in 
this class, systems which display damped oscillations or semi-periodic behaviour. 
One example is oscillations, triggered by DNA damage, in p53, a key protein in- 
volved in cell death and apoptosis pathways [16j El El El • Hormones, such as the 
Human Growth Factor, also show such intermittently periodic behaviour and 
pulsatile secretion [18] . For these systems the recurrent behaviour probably has 
a direct physiological role. In cyanobacteria, for example, various physiological 
processes, like respiration and carbohydrate synthesis, are directly influenced 
by the circadian clock to be in sync with the day-night cycle [TH] . For other sys- 
tems, however, clear oscillatory behaviour is not observed in the wild-type, but 
only in certain mutants. For instance, the NF-kB signalling system, involved 
in immune response in mammalian cells, shows oscillations in nuclear NF-kB 
concentration only in a mutant which contains just one isoform of the NF-kB 
inhibitor, IkBck, and not the other isoforms [3]. Fluorescence measurements of 
the NF-kB system modified so that IkBc* was overexpressed also showed sus- 
tained oscillations over several hours [7]. However, wild- type cells show at best 
damped oscillations 3 . Here it is not clear if the oscillations themselves have 
a significant physiological consequence or are merely a by-product of other re- 
quirements for the wild-type behaviour. Nevertheless, the sub-parts of these 
systems which are potentially oscillatory are important, often essential, compo- 
nents which influence the complex temporally varying wild-type response. 

We will focus on three of the regulatory systems mentioned above: Hesl 
in mammalian embryos, p53-Mdm2 in mammalian cells, and the NF-kB sig- 
nalling system, also in mammalian cells. Fig. [JJ shows the oscillations ob- 
served in experiments for all three. These systems are quite complicated, with 
many components interacting in various ways - including transcriptional activa- 
tion/repression, translation and post-translation regulation, protein-protein in- 
teraction, targeted protein degradation, and active nuclear-cytoplasmic translo- 
cation - composed into a complex network with multiple feedback loops. In this 
review, we mainly describe the mechanisms and structures in these networks 
that allow them to produce the observed oscillations. We will do this using 
simplified mathematical models where the level of description balances the need 
to correctly describe the systems with the need to coarse-grain over some details 
in order to reveal common design features. 

We first begin by defining the overall framework of models we consider, and 
describing the basic ingredients for producing oscillatory behaviour within this 
framework - negative feedback and time delays. In the subsequent sections we 
show how the three examples of Hes, p53 and NF-kB contain these ingredients. 
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Figure 2: Negative feedback loops, (a) A generic multi-species negative feed- 
back loop, (b) the simplest negative feedback loop - a self-repressor, (c) two 
component negative feedback loops for the three examples we examine. In all 
figures a normal arrow represents an activating interaction, and a barred arrow 
represents a repressing interaction. 

Therein we also elaborate, using stability analyses, the requirements for pro- 
ducing oscillations in these systems. Finally, we briefly discuss how to extract 
information about underlying structures from oscillatory time series data. 

Negative feedback 

Fig. [2^ shows a schematic negative feedback loop. The individual nodes in this 
loop are the relevant dynamical variables: they could be protein concentrations, 
gene expression levels, etc. Each variable either activates or represses the next 
one in the loop. By an activator, we mean that if the concentration of protein 1 
goes up it tends to increase the concentration of protein 2 (either by increasing 
its production rate or decreasing its degradation rate). A repressor has the 
opposite effect. Then, a negative feedback loop is simply defined as a loop 
with an odd number of repressors, so that the effect of a perturbation in the 
concentration of any species in the loop eventually feeds back to itself with 
a negative sign. Feedback loops are, in general, the most common network 
motifs in cellular organization, especially when one considers the regulation of 
small molecules [5D]. We concentrate on negative feedback here because it was 
hypothesized by Thomas |21j , and rigorously proven for a wide class of systems 
[22] [23] , that the presence of at least one negative feedback loop is a necessary 
condition for oscillations. 

The simplest negative feedback loop is of course a protein which represses 
itself. There are many examples of such proteins: the main repressor of the SOS 
regulon in E. coli, LexA, also represses its own production [24]; Hesl, mentioned 
above, also represses transcription of its own gene [12] , shown schematically in 
Fig. 0). However, Hesl looks like a one-component loop only if we coarse-grain 
over the intermediate steps in protein production. If we count the Hesl mRNA 
also then it becomes the two component negative feedback loop of Fig. ^p: one 
component, Hes protein, represses the production of the second, Hesl mRNA, 
while the second component activates the first. The other two systems we discuss 
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later in this review have the same loop structure when coarse-grained to the 
two component level. Nuclear NF-kB is known to activate production of IkBq;, 
which inhibits nuclear import of NF-kB by sequestering it in the cytoplasm [3J. 
p53 and Mdm2 work similarly, with p53 activating the mdm2 gene and Mdm2 
sequestering p53 [25] . 

To examine the possible dynamical behaviour produced by a negative feed- 
back loop, we model the dynamics of the concentrations of the components 
using ordinary coupled differential equations. For instance, for the two compo- 
nent loops shown in Fig. [2b: 

dxi 

— = gi{xi,x 2 ) 

dx2 i \ m 

— = g2{xi,x 2 ), (1) 

where x± and xi are the concentrations of the two components X\ and X-i ■ This 
can easily be generalized to longer loops with N components: 

dx ' 

— = giix^Xi-x), i = 1,2, . . . ,N, (2) 

which models a single feedback loop with no cross-links because the rate of 
change of a given variable x$ depends only on itself and the preceding variable, 
Xi-i. In writing such an equation we are assuming that fluctuations in space 
and time are negligible. Thus, there are no stochastic or diffusion terms. The 
functions <?i model both production and degradation of the components, and 
can take many forms depending on the kind of interactions in the system. For 
example, in the p53 example (with X\ — > p53, X2 — > Mdm2), we know that 
p53 binds to an operator site at the promoter for the mdm2 gene and aids 
transcription. Then, under some assumptions the production of Mdm2 can be 
modeled using a term of the form: 

(P/K) h 

l+(p/K) h ' [) 

a sigmoidal monotonically increasing function of the p53 concentration, p (the 
Hill coefficient, h > 1, accounts for cooperativity in the binding of p53 to the 
operator). 

It is difficult to say anything about the behaviour of the general coupled 
differential equation of a negative feedback loop, like Eq. [2j However, it is 
reasonable to constrain the functions gi of Eq. [5] to be monotonic in aff-i- This 
corresponds to saying that a protein that activates a particular process cannot 
change to repress it at some other concentration, and vice versa, which is the 
case for most transcription factors q For monotone systems, not only is there 
no ambiguity about whether the loop implements negative or positive feedback, 



1 In particular, this assumes that the number of proteins bound to the operator site is much 
smaller than the total number of proteins. 

2 Some proteins can both activate as well as repress the same process depending on their 
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but we can also prove rigorously that there is only one fixed point (see Appendix 
A). The question then is whether such a steady state is stable or unstable. For 
two variables, if the fixed point is unstable, and the system's trajectories are 
bounded (again, a reasonable assumption for a biological system) then it must 
show periodic oscillations: this is known as the Poincare-Bendixson theorem 
PH] . For monotone systems this is true even when the loop has more than 
two components |29j , which means that the question of whether oscillations are 
possible boils down to whether the fixed point is stable or unstable. 

Time delays 

Physically, what is required for instability of the fixed point, and hence oscil- 
lations, is a time delay, or a slowing down of the signal going round the loop. 
By signal we mean perturbations of the concentrations away from the steady 
state. If a perturbation in the concentration of one variable instantaneously 
affects the concentration of the next one, and so on, then for a negative feed- 
back loop, any perturbation will be immediately cancelled and the steady state 
will be stable. A sufficiently large time delay on the other hand will produce 
oscillations. This can be readily understood with an example: consider a person 
who walks along a straight line to a given point, marked on the ground. If this 
person is able to take instantaneous decisions, he will approach the mark and 
then stop. This is a stationary solution to the walk kinetics. If, on the other 
hand, it takes some time to realise that the mark has been reached, the person 
will not stop at the mark, but cross it. When eventually the information that 
the mark has been crossed is processed, the person will turn back and walk in 
the opposite direction. The mark will again be reached and overshot, and so 
on. The resulting kinetics will be damped or sustained oscillations about the 
mark. In cellular systems many processes could produce time delays. Table 1 
lists some timescales associated with such processes. 

With the above insight, one can see that a simple way to model oscillations 
using the framework of Eq. [T] is to introduce an explicit time delay into the 
equations: 

gi{xi(t),x 2 (t),xi{t - Ti),x 2 (t- T 2 )) 

g 2 (x 1 (t),x 2 (t),x 1 (t - T 3 ),x 2 (t- r 4 )), (4) 

Of course, this is the most general form, and often it is possible to have fewer 
than four delays. Oscillations are observed even in the simplest case of a linear 
differential equation with a single (sufficiently large) delay (see Appendix C): 

concentration. For instance, CI in lambda phage activates the Prm promoter at low con- 
centrations, but represses it at high concentrations [26]. Another example is the galactose 
regulator GalR, which at high concentrations of galactose activates the promoter galP2 but 
in the absence of galactose forms a DNA loop, which completely represses galP2 1271 . Such 
examples are, however, somewhat rare and we will not consider them. 



dxi(t) 

dt 
dx 2 (t) 

dt 



G 



dx/dt = —x(t — t), which models the 1-component Hes loop of Fig. (2b where 
Hesl represses itself after a time delay r. Although a linear delay rate equation 
is an oversimplification of Hesl production, the physics which lies behind more 
general delay differential equations like Eq. (QJ is the same. The added compli- 
cation is that the functions g\ and (72 are in general highly nonlinear, resulting 
in an amplification of the effect of the delay. 

Putting an explicit time delay like this, of course, does not really shed light 
on the mechanism producing the time delay. There are several possibilities 
which can be used to produce oscillations in a negative feedback loop: 

1. a process that takes a finite minimum time 

2. many intermediate steps 

3. a sharp response by some of the variables 

4. saturated degradation 

5. autocatalysis 

To elaborate: 

1. Rate equations, like Eq. [21 typically model processes which occur with a 
given average rate, such as the binding of a protein to an operator site. A 
hidden assumption is that the time interval between two binding events 
is Poisson distributed, which means that often there is a reasonable prob- 
ability for two events to be separated by a very short time interval (say, 
much shorter than the average time between events). Sometimes, how- 
ever, molecular processes take a certain minimum time. For instance, if 
transcription and translation take a time r after a polymerase binds to the 
promoter, then the rate of production of the protein is more appropriately 
modeled as dx/dt oc P(t — r), where P(t) is 1 if the polymerase is bound 
to the promoter at time t, and zero otherwise. Such logic has been used 
to justify time-delay models in a variety of systems [3TJ [3D1 [32] . This is 
the approach we will take to model the p53 and Hes systems, discussed in 
the subsequent sections. 

2. Processes like transcription and translation have this character because 
they are in fact composed of a large number of intermediate processes: 
the polymerase binds, first forming a closed complex, which then makes 
a transition to an open complex, and then to an elongating complex, fol- 
lowed by many "steps" along the DNA until the polymerase reaches the 
end of the gene. Even if each of these individual steps is a Poisson process, 
the net effect adds up to a time delay. Thus, instead of putting in an ex- 
plicit time delay as in Eq. [U one could work with a negative feedback loop 



7 



with many components. One simple example of such an oscillator is the 
repressilator [33j . which is a negative feedback loop with six components. 

3. The repressilator also has another necessary ingredient - a nonlinearity in 
the gi functions which allows some of the variables to respond to changes 
in the preceding variable in a faster-than-linear fashion. More precisely, 
the repressilator uses sigmoidal functions, like the function [31 to describe 
transcriptional repression, and needs Hill coefficients of at least 2 in order 
to achieve oscillations. The earliest example of a negative feedback oscil- 
lator which uses a nonlinearity like this to produce a sharp response is a 
model by Goodwin [34] where the Hill coefficient needs to be more than 
8. 

4. Such high Hill coefficients are unlikely for biological systems, however. In 
order to get around this "problem" , Bliss, Painter and Marr 35 intro- 
duced another way of producing an effective time delay. They used the 
saturated degradation of one of the concentrations. Saturated degradation 
means that there is an upper limit to the degradation rate of one species, 
thereby allowing it to remain abundant for a longer time, thus effectively 
slowing down the signal travelling around the loop. Such saturated degra- 
dation is quite common in biological systems, especially when proteins are 
tagged for targeted degradation by another protein, as we will show in the 
NF-kB case discussed below. 

5. Finally, autocatalysis, where a molecule activates its own production can 
be used to produce oscillations in systems like Eq. [1] [36] . In fact for two 
variable systems where there is no explicit time-delay this is a necessary 
condition for oscillations [37J . Note that this modification typically makes 
the system non-monotonic. 

This survey of the theoretical requirements for producing oscillations in neg- 
ative feedback loops already allows us to make an interesting observation: Mono- 
tone 2-component loops without an explicit time-delay cannot oscillate, what- 
ever the nonlinearity in the gi functions. We prove this explicitly in Appendix 
A. Thus, if one insists on modelling an oscillating system using two variables, 
one must choose between introducing a time delay, and sacrificing monotonicity. 

2 Ultradian oscillations in biological systems 

We now turn to three biological systems to illustrate these ideas in action. In 
the following, we will briefly give a description of the systems, p53-Mdm2, Hesl 
and NF-kB, in which "ultradian" oscillations have been observed, which have 
time periods of the order of hours, as opposed to "circadian" 24 hour rhythms. 
We will discuss specific details of each system, at the same time emphasising 
that the basic physical mechanism which produces oscillations in all three is the 
same - negative feedback along with time delays. 
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translocation through nuclear pores 38J 


10" 4 s 


diffusion in eucaryotic cell 


1 s 


translation [39] 


30 s 


transcription [39] 


3 min 


mRNA degradation [39] 


3 min 


protein degradation 


10 min to 10 h 


cellular signals [T6][T2]l3] 


lh 



Table 1: Time scales of some cellular processes associated with a single molecule. 
The upper part of the table indicate the processes which are usually neglected in 
writing the rate equations of regulation newtworks, while the lower part indicate 
proceeses which are usually accounted for. 




Figure 3: A sketch of the feedback loops controlling the concentration of p53 
(a), Hesl (b) and NF-kB (c). 
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2.1 P 53-Mdm2 



The protein p53 is responsible for inducing apoptosis in cells with damaged 
DNA [5D]. The concentration of p53 is usually kept low by a feedback mech- 
anism involving another protein, Mdm2, which binds to p53 and promotes its 
degradation. When the DNA is damaged, the cell expresses a number of kinases 
which phosphorylatc SER20 in p53, changing its affinity to Mdm2. This results 
in oscillations in the concentration of p53, observed both in western blot analy- 
sis [THIHI] (cf. Fig. [Ha)) and in single cell fluorescence experiments [HIS]- The 
standard explanation for the overall increase in the concentration of p53 is that 
its phoshoprylation decreases its affinity to Mdm2, shifting the thermodynamic 
equilibrium towards higher concentrations. 

Apart from not explaining the oscillations, this argument does not agree 
with other experimental evidence. Equilibrium isothermal titration calorime- 
try experiments have shown [42] that phosphorylation at SER20 decreases, 
not increases, the dissociation constant between p53 and Mdm2 from kr> = 
575 ± 19 nM to krj = 360 ± 3 nM. The same effect is observed in vivo [33], 
where p53ASP20 (a mutated form which mimics phosphorylated p53) binds 
Mdm2 more tightly than p53ALA20 (which mimics unphosphorylated p53). 
Moreover, single cell experiments [5] show a slight decrease in the concentration 
of p53 after DNA damage, which cannot be explained by the standard argument. 

We showed that a simple model of the p53-Mdm2 system that incorporates 
the time delay associated with some relatively slow processes within the cell can 
account for the experimental facts in a simple way [30] . The feedback mechanism 
is sketched in Fig. C^a) and the associated time-delayed rate equations are 

S — a ■ pm — b ■ p (5) 

p(t — t) — pm(t — t) 
c- — d ■ m 

kg + p(t — t) — pm[t — t) 

- (^(p + to + k) — y / (p + m + k) 2 — Ap ■ rnj , 

where the delay r takes into account the half-life of mRNA, the diffusion time, 
the time needed to cross the nuclear membrane and the transcription/translation 
time. One can solve Eqs. ([5]) numerically, simulating the damage to DNA 
by a sudden change in the dissociation constant k. Fig. [4] shows the onset 
of oscillations in this model in response to a sudden decrease in k at time 
t = 2000 s. 

Several observations emerge from the simulation results in Fig. [4[ Most 
interestingly, what triggers oscillations is a decrease in the dissociation constant 
fc, while any increase in its value just shifts the equilibrium towards higher values 
of p (cf. inset of Fig. [4]). Moreover, just after t — 2000 s the concentration of 
p53 decreases, as observed in the experiments. Finally, the concentration, pm, 
of the complex p53-Mdm2 is, at any time, essentially identical to the minimum 
of p and to (gray curve in Fig. [3]), indicating that the complex is saturated. 



dp 
dt 
dm 

~dt 

pm = 
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Figure 4: Oscillations displayed by the numerical solution of the dynamic equa- 
tions of p53. The numerical values of the rates are a = 3T0~ 2 s _1 , b — 10~ 4 s" 1 , 



1 s- 



1 s- 



S = 1 s- 



180, k 



28 and r = 1200 s. At time 



c = 

t = 20000 s the value of k is decreased of a factor 10. In the inset, the solution of 
the same equations, where the value of k is increased of a factor 10 at t = 20000 
s. The equations are solved with the Adams algorithm. 
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Figure 5: The height of the response peak Ap with respect to the quantity that 
multiplies k, mimicking the stress. The dotted line indicates that the system 
does not display oscillatory behaviour. 

The relative change, Ap, in the maximum concentration of p53 reached after 
the simulated stress event is displayed in Fig. [5] as a function of the change in k. 
The value of Ap displays a sigmoidal behaviour if the stress decreases the value 
of k. In contrast, in the absence of a time delay, this curve is linear with respect 
to k |30j . Thus, the time delay and the oscillations are crucial for producing a 
sensitive system that can respond to stress in a sharp, faster-than-linear manner. 

We also undertook a detailed analysis of the response of the oscillations to 
changes in the parameters [30] . It turns out that the oscillations do not change 
qualitatively when parameters a, b and c are varied (by upto five orders of magni- 
tude), whereas a decrease in d or k g suppresses the oscillating behaviour. These 
parameters are associated with, respectively, the the degradation of Mdm2 and 
the binding of p53 to the DNA. One could speculate therefore that these pro- 
cesses play an important role in the production of tumors that arise when the 
p53 response mechanism fails. This speculation is, in fact, consistent with the 
observation that around the 45% of all tumors display mutations in the p53 
region which binds to DNA [4"3] . 

2.2 Stability analysis of the p53-Mdm2 model 

The overall behaviour of the proteins with respect to time evidently depends on 
the parameters of the dynamic equations, that is the production and degradation 
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rates and the delay. Over long timescales, the protein concentrations can either 
converge to a steady state or a limit cycle (i.e., sustained oscillations). Chaotic 
behaviour is never observed within this model. 

A careful analysis of how the dynamics of the protein concentration depends 
on the parameters of the rate equations has been done by Neamtu and coworkers 
in ref . [35] . The main conclusion of this work is that under the condition that the 
dissociation constant k between p53 and Mdm2 is small, there exists a critical 
delay, tq, above which the system shows oscillations (see details in Appendix 
C). 

In the language of dynamical systems, the transition when the delay, t, 
crosses r is a Hopf bifurcation. A Hopf bifurcation is the generic type of 
bifurcation that occurs when a stable fixed point of the system becomes unstable 
and turns into a limit cycle, a dynamically oscillating state, as consequence of 
a change in some parameter of the system (see Appendix B for more details). 

2.3 Hesl and its mRNA 

The transcription factor Hesl controls the differentiation of neurons in mam- 
malian embryos [3S]- Its concentration is controlled by a feedback loop built 
out of Hesl and its own mRNA (see Fig. [3]). Like p53, its concentrations dis- 
plays oscillations when the cells are stimulated with serum [12]. The period is 
similar to that of p53, approximately two hours, and the oscillations last for 
« 12 hours. Hesl is known to repress another protein called Mashl. Both the 
knocking out of Mashl and its continuous expression by means of retroviral in- 
troduction results in a lack of cell differentiation. Only when there are periodic 
oscillations in the concentration of Hesl (and thus of Mashl) is there proper 
differentiation of neuronal cells [46] , showing that the temporal variation of the 
protein concentrations are critical for the developement of the nervous system. 

The relation of Hesl oscillations to segmentation and spatial patterns has 
been studied by several authors. It is well known that oscillations in Hesl is 
part of Notch signalling and Hirata et al Ref. [12] studied the coordinated 
somite segmentation in the presomitic mesoderm of mice embryos and found 
a correlation between the oscillations of Hesl and initiation of somites. The 
general issue of segmentation in vertebrates has been further studied in Ref. 
[TBI [HI [T5] indicating that the oscillations in the Notch pathway signaling are 
intimately related to the Wnt signaling. At the caudal end of the embryo Notch 
and Wnt oscillates out of phase when the gradient in the Wnt level is over a 
given threshold |14[ 115] . New cells are provided and move away from the caudal 
end each setting a boundary for a somite segmentation. This continues until 
all somites are produced which subsequently gives rise to the spinal cord |13j . 
Ultradian oscillations of signaling pathways are thus extremely important for 
segmentation. 

The control mechanism of Hesl oscillations again involves a negative feed- 
back loop with one activation and one repression: the transcription of the mRNA 
of Hesl activates the production (translation) of the protein, and Hesl represses 
the transcription of its own mRNA. The main difference compared to p53 is that 
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here one of the nodes represents an mRNA, not a protein. However, this differ- 
ence is merely semantic: the control network still has two nodes, of which one 
is activating and one is repressing. The molecular species of these nodes are 
immaterial to the description of the oscillations 0. 

We model the system using the following equations (cf. [5T] ) 



where s and r are the concentrations of Hesl and its mRNA, respectively. The 
meaning of these equations is that mRNA is produced at a rate a when Hesl is 
bound to the DNA. The probability that Hesl is bound to DNA is k h /(k h + s h ), 
where k is a characteristic concentration for dissociation of Hesl from the DNA, 
and h is the Hill coefficient that takes into account the cooperative character 
of the binding process. k r and k s are the spontaneous degradation rates of the 
two proteins, while r is the delay associated with the molecular processes that 
we do not want to describe explicitly (transcription, translocation, etc.). Ref. 
|12j suggests that r rna and r^esi are of the order of 25 minutes. The value of 
the time delay is difficult to assess, since it is determined by a combination of 
various molecular processes. One can guess that its order of magnitude is tens 
of minutes. 

The numerical solution of Eq. (JT]) is displayed in Fig. [6l The oscillations 
have a period At rs 170 min. The dependence of the time period on the delay r 
is shown in Fig. [7l For any delay in the range 10 < t < 50 min, the oscillation 
period is consistent with that found experimentally, and so is the time difference 
between the peaks in Hesl and mRNA, w 18 min. For r < 10 min, the system 
shows no oscillations. To check the robustness of the results, we have varied a, 
j3 and k over 5 orders of magnitude around the basal values, listed in the caption 
to Fig. [6l and observed no qualitative changes to the oscillatory behaviour. On 
the other hand, a decrease of k s and k r disrupts the oscillations. This is because 
these two quantities set the timescale of the dynamics, Increasing this timescale, 
keeping r constant, is equivalent to decreasing r and thus no oscillations are 
produced. 

The behaviour of Hesl is very similar to that of p53, both in the features of 
the oscillations and in the lag delay before they start. This is not unexpected, 
since the structure of Eq. ([7]) is very similar to the structure of Eq. (JSJ) and to 
any time delay equation describing a two component negative feedback loop. 

3 We could of course introduce two more nodes in the p53 network representing the mRNA 
of p53 and Mdm2. However, this would not change the logic of the loop, replacing an activation 
by two successive activations. As discussed earlier adding more intermediate nodes introduces 
a time delay. Since the p53 equations had an explicit delay it seems redundant to add the 
mRNA nodes also. 
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Figure 6: The concentration of Hesl (solid curve) and uts mRNS (dashed curve) 
obtained solving numerically Eq. |7|). 
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Figure 7: The oscillation period At of Hesl as a function of the delay r. 
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2.4 NF-kB and IkB 

The NF-kB family of proteins is one of the most studied in the last ten years, 
being involved in a variety of cellular processes including immune response, in- 
flammation, and development [31117]. NF-reB can be activated by a number of 
external stimuli [35] including bacteria, viruses and various stresses and pro- 
teins (for instance, refs. [3[ [JJ used the tumor necrosis factor-a, TNF-a). In 
response to these signals it targets over 150 genes including many chemokines, 
immunoreceptors, stress reponse genes, as well as acute phase inflammation re- 
sponse proteins [48] . Each NF-kB has a partner inhibitor called IkB, which 
inactivates NF-kB by sequestering it both in the nucleus as well as in the cy- 
toplasm. In fact, the IkB proteins come in several isoforms a,(3,e [3 [49], and 
perhaps others |47J. Some of these isoforms are, in turn, transcriptionally ac- 
tivated by NF-kB, thus forming a negative feedback loop which is essentially 
identical in structure to the other two discussed above. 

The potential for this negative feedback loop to produce oscillations in 
the nuclcar-cytoplasmic translocation of NF-kB was initially shown by elec- 
trophoretic mobility shift assay experiments [3 . They found that wild-type cells 
and mutants containing only the InBa isoform showed damped oscillations. In 
contrast, cells with only the IkB/3 or e isoforms do not show oscillations. This 
conclusion was bolstered by single-cell flourescence imaging experiments which 
show sustained oscillations of nuclear NF-kB in mammalian cells [JJ, with a 
time period of the order of hours. In these experiments InBa was overexpressed 
hence the system behaves like the mutant which has only the IkBch isoform. 

As sustained oscillations are clearest here, we will focus our modelling on this 
mutant. The following cellular processes, summarized in Fig. [8j are important 
for this system on the timescales we are interested in (i.e., we ignore processes 
which are very slow): 

• NF-kB, when in the nucleus, activates transcription of the InBa gene 
(henceforth we will drop a unless we are explicitly talking about more 
than one isoform), producing IkB mRNA in the cytoplasm. 

• The IkB mRNA is translated to form IkB protein. 

• The IkB protein can be transported in and out of the nucleus. 

• In both compartments, IkB forms a complex with NF-kB. 

• The NF-kB-IkB complex (henceforth referred to as {NI}) cannot be im- 
ported into the nucleus. However, if it forms within the nucleus it can be 
exported out. 

• Free NF-kB behaves in exactly the opposite way. Free NF-kB is actively 
transported into the nucleus but not from the nucleus to the cytoplasm. 

• The cytoplasmic {NI} complex is tagged by another protein, the IkB ki- 
nase (IKK), for proteolytic degradation. This results in degradation of 
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Figure 8: The important interactions in the NF-kB system. Green arrows in- 
dicate transcription and translation. Blue arrows indicate transport processes. 
Purple double arrows indicate complex formation. The red dashed arrow indi- 
cates IKK triggered degradation of IkB when complexed to NF-kB. 
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IkB only, releasing NF-kB. Note that this degradation does not occur for 
free IkB. 



On the timescales of interest, there is no net production or degradation of 
NF-kB. It simply cycles in and out of the nucleus, i.e., the sum of nuclear and 
cytoplasmic NF-kB concentrations is a constant. Amongst the above listed pro- 
cesses, the association and dissociation of the complex {NI} occurs fast enough 
that the concentration of the complex can be taken to be always in equilibrium 
with the free NF-kB and IkB concentrations. This allows us to describe the 
system using a very simple model consisting of only three variables [50] , nuclear 
NF-kB (N n ), cytoplasmic IkB (J) and IkB mRNA (I m ): 

dNn = A (^N A _ B IN^ 
dt e + I 5 + N n ' K ' 

d±_ c (l-N n )I 
dt ~ Im C e + I ■ {W) 

The terms in equation (JSJ model the nuclear import and export of NF-kB. 
Eq. ([9]) models NF-KB-activated transcription of the IkBq: gene and sponta- 
neous degradation of the mRNA. Finally, eq. (|10|) has terms for translation of 
IkB mRNA into the protein and IKK mediated degradation of IkB. The exter- 
nal signal is supplied by IKK that enters the equations through the parameter, 
C, which is proportional to IKK concentration. 

These equations produce sustained oscillations in all variables over a wide 
range of parameter values. Fig. [5] shows the result of simulations which use 
parameter values from ref. [3J. The oscillations produced by this model match 
the experimentally observed features well, in particular, the shape, phase, time 
period and frequency response are correctly reproduced [50] . 

Note that this is certainly not the only model able to reproduce the experi- 
mental oscillations observed in NF-kB. Hoffman et al. have constructed a long 
list of chemical reactions between 26 different molecules in the NF-kB system, 
including 65 numerical parameters O [49] . The three variable model above was 
formed by reducing this larger model [SOj, also producing a 7- variable and 4- 
variable model in the process. Hayot and Jayaprakash have also built a model 
with seven variables for NF-kB oscillations |51j . Another model, which is sim- 
ilar to Hoffman et al's model, but has an additional feedback loop is described 



2.5 Saturated degradation of IkB 

A key element in the model is the saturated degradation of cytoplasmic IkB in 
the presence of IKK (second term in Eq. (fTUj) V This saturation occurs because 
the level of the NF-kB-IkB complex saturates, and this complex is needed for 
IKK triggered degradation of IkB. A stability analysis of the system shows the 
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time(min) 

Figure 9: Oscillations of nuclear NF-kB (N n ), red, and cytoplasmic IkB, green, 
for A = 0.007, B = 954.5, C = 0.035, S = 0.029 and e = 2 x 10" 5 (these 
parameter values are derived from the ones used in ref. [3], see |50j.) In order 
to facilitate comparison with the experimental plot, the x-axis has been limited 
to 600 minutes, but the oscillations are in fact sustained (see Fig. [TTk). 
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importance of the saturated degradation for oscillations. We begin by examining 
the fixed points of the system. 

The fixed point values of N n , I m and / are solutions to 



A 



(1 - Nn) 
e + I 



B- 



IN n 
S + N n 



0, 



N 



0. 



(1 - N n )I 

Im - C± ^- = 0. 



I m and I can be eliminated using I m = N 2 , giving I — (N 2 e)/(C — CN n — 
N 2 ). From this we find that the fixed point value of N n is a solution of the 
equation (C - CN n - N 2 ) 2 = {BCe 2 N*)/[A(6 + N n )], or equivalently, 

B 



N>+(5+2C)N*+C 



N*+C[(C-2)6-2C}N 2 +C 2 (l-2S)N n +C 2 S = 0. 



In general, this has two real solutions, one with C — CN n — N 2 > and the 
other with C — CN n — N 2 < 0. The latter results in a negative value for I and 
therefore is not an acceptable solution. Thus we are left with only one fixed 
point. 

Next, we linearize the equations around the fixed point, which gives the 
Jacobian (see also Appendix B) 



V 



SBI 

2N n -1 



- 



ATI 



CI 
e+I 



1 





Ce(l- 



BN n 



This matrix can be used to examine the stability of the fixed point: it is 
unstable if any of the eigenvalues have a positive real part. For the NF-kB 
system, Fig. [10] shows how the stability depends on the parameter e. When e 
is small compared to the steady state value of J, the degradation rate of IkB is 
independent of /, which is what we mean by the term "saturated degradation" . 
On the other hand, when e is large then the degradation is proportional to I 
and is not saturated. Fig. [TU] shows that, indeed, the fixed point is unstable 
when e is small compared to I. Physically, this is because saturated degrada- 
tion introduces an effective time-delay into the feedback-loop: it allows IkB to 
accumulate and stay around longer than with non-saturated degradation. 



2.6 Spikiness and sharp response of the NF-kB oscillation 

One property of the oscillations of nuclear NF-kB, in figured! that stands out is 
that they are extremely spiky. By choosing different parameter values one can 
also get soft oscillations (see Fig. [TTk.b). however, for biologically relevant 
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Figure 10: (A) The stationary value of J as a function of the parameter e which 
controls the binding between NF-reB and IKK in the cytoplasm. (B) and (C) 
The projection of two trajectories at different values of e onto the N n — I plane. 



parameters the oscillations are spiky. We suggest the following measure of 
spikiness: 

z= max-min 
average 

In other words, Z is the ratio of the amplitude of the oscillations to the average 
level. Then, we will call a particular oscillation spiky if Z > 2 and soft otherwise. 
Further, Z = indicates the absence of oscillations. Fig. [TTh shows a contour 
plot of Z values on the B — C plane. In some directions the system transitions 
quickly from no oscillations to spiky oscillations, wheras in other directions there 
is a softer transition. In general, the existence and spikiness of the oscillations 
is very robust to changes in most of the parameters of the model [SO] . 

There is one parameter, however, to changes in which the system shows a 
very sensitive response. That parameter is the external input, IKK. Fig. [12] 
shows that both the spike height (or peak level), as well as the spike duration, 
can change by large amounts in response to small changes in IKK level. Notice 
that this sensitivity is particularly high in IKK ranges which are near the transi- 
tion from spiky to soft oscillations. One way of quantifying this sensitivity is to 
measure the expression level of a gene whose transcription is activated by NF- 
kB. Imagine a gene whose upstream regulatory region contains a binding site for 
NF-kB dimers. Upon, binding the gene promoter is activated: G+2N ^ k ° o n G* . 
Experimental measurements of NF-«:B-dependent gene expression suggest that 
many genes closely follow the oscillations of NF-kB [5]. This corresponds to 
the case where the binding of NF-kB to the operator is in equilibrium, i.e., k on 
and k ff are much larger than the rates of all other processes in the NF-kB 
system. In that case the gene activity, G* , will follow the NF-kB concentration: 
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Figure 11: (A) Spiky oscillations in the NF-kB model (parameter values are 
identical to those used to produce Fig. [9]). (B) Soft oscillations, produced in the 
same model using a larger value of e, keeping all other parameters unchanged. 
(C) A countour plot of the Z measure of spikiness (see text), in the B — C 
parameter plane (other parameters are unchanged). The Z = 0.01 contour 
maps, approximately, the region of oscillations. The Z ~ 2 contour shows the 
region of spiky oscillations. The black dot marks the value of B and C used in 
(A) and (B), as well as in Fig. [9l 
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Figure 12: Sensitivity to IKK. A. Spike duration, the fraction of time N n spends 
above its mean value, as a function of IKK concentration. B. Spike peak, the 
maximum concentration of nuclear NF-kB, as a function of IKK concentration. 
In both plots, the black dot shows the IKK value used in Fig. [51 while blue and 
red signify, respectively, regions of spiky and soft oscillations. 



G* = k /fc" + jy2 • In this case the downstream gene inherits the high sensitiv- 
ity of the NF-kB signal to IKK: the peak gene activity as a function of IKK is a 
very steep sigmoidal curve with an effective Hill coefficient of over 20 [50] . Such 
a large value is very unusual for biological systems, and is much larger than the 
values obtained by other mechanisms |53[ 154] . 



3 A tool to analyse oscillation patterns in time 
series 

After having discussed the mechanism underlying the production of oscillations 
and studied in detail three examples of "ultradian" oscillations in cells, we would 
like to describe an easy tool [55 capable of providing information on the archi- 
tecture of the underlying network, given the experimental time series data for 
oscillating biological systems. 

More precisely, given the time sequence of maxima and minima of the con- 
centrations of the species during the oscillations, the method allows one to assess 
whether the observed sequence is compatible with a a single underlying negative 
feedback loop, describable by Eq. O If it is, the method further predicts the 
logical structure of the loop, i.e. the order of the proteins within the loop, as 
well as which protein acts as an activator and which as a repressor. 

The method is grounded in a mathematical result valid for all monotone 
systems describable by Eq. O For such a system, we can prove the following 
statements about the sequence of maxima and minima in the time series: 
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Figure 13: Reconstructing the underlying loop from the time series. The time 
series shows p53-Mdm2 oscillations observed in single cell fluorescence experi- 
ments [5]. Using the rules mentioned in the text, the sequence of maxima and 
minima is converted to the feedback loop shown at bottom right. 

• the extrema (maxima or minima) have to follow the order of the variables; 
for example, after a maximum of variable i — 1 there can either be a 
maximum or a minumum of variable i, and 

• if two successive extrema are of the same kind (both maxima, or both 
minima), then variable i is activated by variable i — 1. Conversely, if they 
are of opposite nature (one maximum and one minimum), then variable i 
is repressed by variable i — 1. 

These statements can be used to reconstruct the structure of the underlying 
loop from the time series, if it is compatible with the above restrictions. We 
illustrate this using the example of p53-Mdm2 time series shown in Fig. [T3l 
The reconstructed loop is exactly the structure we used in Eq. [5j p53 activates 
Mdm2, which represses p53. Of course, in this case, there are only two possible 
loop structures, and the correct interactions are already known experimentally 
so the information gained is minimal. The method comes into its own for sys- 
tems with more variables. One example is the cyanobacterial circadian clock. 
Fig. HH shows oscillations in the expression of three circadian clock genes in 
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Figure 14: Reconstruction algorithm applied to circadian rhythms in cyanobac- 
teria. Data from [lOj . The algorithm predicts that kaiA activates kaiCl, which 
represses kaiBS, which, in turn, activates kaiA. 

Synechocystis and the reconstructed loop. Two of the interactions have been 
experimentally observed in another strain of cyanobacteria [56j , but the activa- 
tion of kaiA by kaiBS is a new prediction. 

Notice that there are many possible maxima/minima sequences which are 
not compatible with the above rules. For example, in order to be compatible 
with a single loop, each species must have exactly one maximum and one min- 
imum during one period. Moreover, if the method predicts an even number of 
repressors, one should consider an oscillating positive feedback loop, which is 
impossible (see [55] for a full list of non-allowed cases). In all these situations, 
one should conclude that the mechanism causing the oscillations is more com- 
plicated and cannot be reduced to a single feedback loop. Even if the structure 
of the real protein interaction network is more complicated that a single loop 
(which is usually the case), the method can hint at the interactions which are 
most relevant for the oscillating mechanism and help in building up a zero-order 
model of the system. 

Finally, we note that while the method is mathematically rigorous for sys- 
tems without time delays, it works even if there are unobserved chemical species 
taking part to the loop 55J. As intermediate species introduce time delays it 
is likely that the result can be extended to systems with an explicit time delay 
added to Eq. H 
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4 Summary and Outlook 



Questions about cellular processes that use complex, temporally varying signals 
can be crudely classified into the 'How' and 'Why' groups. 'How' questions 
deal with the structure of the regulatory network and the range of dynamical 
behaviour it can produce: how does the network produce oscillations? what 
are the necessary and sufficient mechanisms? how are they implemented in 
real cellular systems? 'Why' questions deal with the physiological role of the 
particular dynamical behaviour produced by the network: why does p53 start 
oscillating in response to DNA damage? do oscillations carry some information 
that can be decoded by downstream genes? could a non-oscillating system have 
worked equally well for the cell? 

In this review we have mainly investigated the 'How' questions for one sub- 
set of cellular response systems that use temporally varying signals: eukaryotic 
systems exhibiting "ultradian" oscillations. The three systems we have mod- 
elled - p53, which is important for apoptosis, Hesl, which is part of the Notch 
cycle responsible for somite segmentation, and NF-kB, a key protein in immune 
response - all turn out to have the same basic design that produces oscillations. 
The two key aspects of this design are the presence of negative feedback and time 
delays. There are many ways of producing an effective time delay in cellular 
processes. In particular, we emphasised, using the example of NF-kB, saturated 
degradation as one such mechanism. 

Saturated degradation plays a role in many models of oscillatory negative 
feedback-loops. The earliest use of this mechanism is in the model of Bliss, 
Painter and Marr [35] , which is similar to our NF-kB model. Saturated degra- 
dation has also been used by Goldbeter in various models of cellular oscillations, 
e.g., the cell cycle [57], development in myxobacteria [53], yeast stress response 
[55] and the mammalian circadian clock [50] . It is also found in models of cal- 
cium oscillations in cells [5TJ [55]. In the p53-Mdm2 system too, as described 
earlier, there is a similar saturated complex formation, where one component 
of the complex is subsequently degraded. This is so similar to the NF-kB case 
that it should be possible to construct a model for p53 oscillations which has 
ordinary differential equations without an explicit time delay (see for instance 
[5J). Because saturated degradation is easily implemented by complex formation 
we conclude that it would be a very useful general mechanism for producing an 
effective time delay in cellular systems. 

Although we have focussed more on how the oscillations are produced, the 
'Why' question of whether the oscillations have a direct physiological role in 
these systems is of course an important one. In NF-kB, oscillations are observed, 
not in wild-type cells, but rather in mutants or cells which overexpress certain 
proteins. Therefore, it is quite possible that oscillations themselves play no 
direct role in the wild- type response [63] , However, wild-type cells do show a 
complex temporal variation of NF-kB of which the spiky response of the luBa 
module is an extremely important component (luBalpha is the only isoform 
whose knockout mutants are not viable [3]). Hesl oscillations, even though 
damped, might have a direct relevance for spatial pattern formation in embryos. 
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Hesl is a key element in the Notch oscillating cycle, which, in mouse and chicken 
embryos, is coupled, out of phase, with the Wnt cycle. In each period of these 
cycles, a somite is formed, thus leading to vertebrate segmentation. Further 
work should clarify whether the Hesl oscillations drives these cycles, or is slaved 
to them, or has a different role entirely. In the case of p53, the oscillations 
appear to be sustained, therefore it is likely they are directly important for the 
apoptosis pathway, but again it is not clear precisely how. 

Even in the absence of unambiguous evidence that oscillations have a direct 
functional importance, investigations into the mechanisms underlying the oscil- 
lations have implications for some of the 'Why' questions. For instance, noting 
that an oscillating signal carries much more information than a steady signal, 
[17l [64] suggested that differential regulation of downstream gene circuits could 
be achieved if they were sensitive to the frequency of the oscillation. Indeed, as 
we have shown [50], it is easy for a single downstream gene to act as a 'low-pass 
filter'. Regulatory circuits with multiple genes could be constructed to have 
particular frequency responses. However, it is perhaps more fruitful to look for 
a physiological role in other properties of the dynamics, rather than the peri- 
odicity. One such property, especially prominent in NF-kB oscillations, is the 
spikincss. A spiky signal of this kind carries even more information than a soft 
oscillation or a steady state level. Spiky pulses, whether periodic or random 
or isolated, are in a sense less "expensive" : If a downstream gene is expressed 
when p53 crosses some threshold it is energetically and metabolically cheaper 
to have a spike whose peak crosses the threshold for long enough to trigger the 
gene, than to produce and maintain the protein above that level for a long time. 

In this context it is noteworthy that several of the most ubiquitous signalling 
molecules in eukaryotes, hormones, exhibit recurrent spikes. The spikes some- 
time come in a nearly periodic fashion but may also appear more randomly 
Why hormones appear in spikes is also not known. Again, we speculate that 
spikes are an efficient way to trigger other regulatory mechanisms in the body 
as compared to a steady level which do not deliver 'kicks' to stimulate other 
compartments. An added benefit is that the time between spikes allows an extra 
level of regulation. In addition, it might not be good for a biological organism 
to be subjected to a constant level of hormones all around the clock. 

A related issue where spikiness has implications is that of cross-talk between 
signals. Thus, a high constant level of a protein may potentially interfere with 
other signals being sent, for instance by saturating common receptors. In con- 
trast, a 'quantization' of the signal into spikes allows different signals to use 
common receptors without fear of interference. An analogy to cars is perhaps 
useful to visualize this difference: If there were no traffic lights, it would still be 
easy to drive a car to your destination provided the traffic was not too heavy. 
However, if all cars were made 10 times longer it would become much harder to 
drive at the same density of traffic. Indeed, railway tracks have less intersections 
because trains are, in effect, very long cars. 

On a more concrete level, our analysis of the NF-kB network suggests that 
the spikiness of the oscillations was correlated to a sharp, threshold-like response 
of the system. The output (nuclear NF-kB concentration), we found, could be 
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extremely sensitive to changes in the input (IKK) while remaining relatively 
insensitive to other parameters; a clear prediction that awaits experimental 
verification. Such a threshold-like response is likely to be very important for 
the differential regulation of downstream genes. It would be very interesting to 
know if there is a causal connection between spikiness and sharp responses in 
this system, as well as other regulatory systems. 

Overall, we conclude that oscillations, especially spiky ones, have many use- 
ful properties that could be exploited to improve the speed and efficiency of 
signalling and response systems. 

We are grateful to Alexander Hoffmann, Terry Hwa, Arnie Levine, Eric 
Siggia, Galit Lahav and Ian Dodd for interesting discussions on oscillating ge- 
netics. S. K., M. H. J and K. S. acknowledge support from The Danish National 
Research Foundation and Villum Kann Rasmussen Foundation. G. T. acknowl- 
edge support from the FIRB 2003 program of the Italian Ministry for University 
and Scientific Research. 

A Properties of monotone systems 

In this section we study the fixed point properties of a feedback loop composed 
of an arbitrary number, N, of nodes whose dynamics is given by Eq. ^ in the 
main text, which we repeat here: 

< ^=9 < f' R \x i ,x i - 1 ) i = l...N. (12) 

Our analysis proceeds by noting that, using the monotonicity condition, 
we can write explicit functional relations between neighboring variables in the 
steady state (when dxi/dt = 0): 

gl A ' R \x*,xU)=0 =* x! = ti A>R) (xU) (13) 

Notice that the functions fi have the same monotonicity properties as the <^s 
with respect to the second argument (for this it is necessary that gi(x,y) be a 
monotonically decreasing function of x). By iterative substitution, we obtain: 

x * = /iOi-l) = fi(fi-l(%i-2)) = ■■■ = 

= fi o fi-i o/ w o...o f i+ i(x*) = Fi(x*) (14) 

where o denotes convolution of functions. Here, we introduced the function 
Fi{x), which quantifies how the species i interacts with itself by transmitting 
signals along the loop. Notice also that if Eq. (fT4"|) holds for one value of i, then 
it holds for any i, since it is sufficient to apply fi+i() on both sides to obtain 
the equation for x* +1 and so on. For feedback loops, much useful information 
can be obtained from the properties of Fi(x). Firstly, by applying the chain 
rule, we obtain the slope of Fi{x) at x: F((x) = Y\j fj{xj)\ Xi=x . The r.h.s is 
always greater (less) than zero if the number of repressors present in the loop is 
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even (odd). In the former case, there can be multiple fixed points, i.e., this is 
a necessary condition for multistability. On the other hand, when there are an 
odd number of repressors, then Fi(x) is positive and monotonically decreasing, 
meaning that there is one and only one solution to the fixed point equation 

x* i =F i (x* i ). 

To perform the stability analysis, we write the characteristic polynomial 
evaluated at the fixed point: 

Y[ [A - d x gi(x,y)\ x = x ,] = \\_dyg t {x,y)\ x=x * . (15) 

i i 

The above equation can be greatly simplified using the relation F'(x) = 
Y\i d y gi{x, y)/d x gi(x, y), which is a consequence of the implicit function theorem 
and the chain rule. One then obtains the following equation: 

n no as) 

where the hi = —d x gi(xi, Xi = i)\ x * are the degradation rates at the fixed point. 
Notice that, because F'(x) is always negative in a negative feedback loop, all co- 
efficients of the characteristic polynomial are non- negative, hence it can not have 
real positive roots. This means that the destabilization of the fixed point can 
only occur via a Hopf bifurcation, i.e. with two complex conjugate eigenvalues 
crossing into the positive real half-plane. 

In the simple case in which all the degradation rates are equal and unchang- 
ing (i.e, hi = 7, a constant) the roots of the polynomial (|16[) in the complex 
plane are the vertices of a polygon centered on —7 with a radius \F'\ as sketched 
in Fig. [15] Therefore, the fixed point will remain stable as long as 

|F'(a;*)|cos(7r/JV) < 7. (17) 

In this case, Hopf 's theorem (see Appendix |B|) ensures the existence of a periodic 
orbit close to the transition value, whose period is: 

T = 2n/Im(X) (18) 

which, in the simple case of equal degradation timescales, becomes T — 2ir/[\F'(x*)\- 
sin(7r/iV)]. Notice that the Hopf theorem does not ensure that the orbit is sta- 
ble; however, since the system is bounded and there are no other fixed points, 
we expect the orbit to be attracting, at least close to the transition point. 

Notice that condition (fP7|) is always satisfied when N = 2. This result also 
extends to the more general case where degradation rates are unequal. Thus, we 
have proven that a two-component monotone negative feedback loop without 
an explicit time delay can never show oscillations. 

B Hopf bifurcations 

Usually, the bifurcation takes place under the variation of some external 'control' 
parameter \x and the new state appears at a critical value fi c of this parameter. 
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Figure 15: Sketch of the Hopf bifurcation in the eigenvalue complex plane, in 
the case in which all the degradation rates are equal to a constant 7. 

For genetically regulated networs, the control parameter could be some external 
chemical stimuli, a production rate, or a binding constant, just to mention a few 
examples. To illustrate a Hopf bifurcation let us first consider a completely gen- 
eral two-dimentional dynamical system (without delay) defined by two coupled 
nondinear differential equations 

( m ) ~ ( g(p,m) ) ^ 

(we use variables (p, m) to resemble p53 and mdm2). The stationary fixed point 
(p* 7 m*) is determined by 

f(p*,m*)=g(p*,m*)=Q (20) 

Using standard routines one linearizes around the fixed point by means of the 
Jacobian matrix 

/ df df_ \ 
J = t X ) (21) 

\ dp 9m /, 

where the star symbolizes that we insert the fixed point into the Jacobian. The 
resulting set of eigenvalues Ai,A 2 can either both be real or be a set of two 
complex conjugates. This is the case we are interested in here: 

Xi = a + iu> , A 2 = a — iu; (22) 

For a < (correpsonding to /1 < ji c ) the fixed point (p*,m*) is stable whereas 
the limit cycle becomes stable for a > (corresponding to \i > n c ) and per- 
forming oscillations on the limit cycle with a frequency u> defining th period 
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of the oscillation. For a genetic system this period could for instance be cir- 
cadian, ultradiand or related to cell cycles. Just above the bifurcation, the 
limit cycle can be well approximated by a circle whoes radius for super-critical 
Hopf bifurcations generally grows continuosly as y ' \x — fj, c (as in contrast to a 
sub-critical Hopf bifurcation where the radius will jump and exhibit hysteresis 
effects). When the control parameter \i is much larger than [i c the limit cycle 
may easily be deformed in various ways 'away' from the circle. 

When we introduce a time delay in the equations, as discussed above, we 
can formally write it by introducing the delayed variable p T into the equations: 



V \ = ( f(p,m) 
rh J \ g(p,m,p T ) 



(23) 



We find numerically that the system still exhibits Hopf bifurcations when cross- 
ing critical lines in the parameter space. In fact, we find that one can drive the 
system through Hopf bifurcations by increasing the value of the time delay t 
through specific values. 

The fixed point ofp is of course equal to the fixed point ofp T , i..e.: f(p*, to*) = 
g(p* , m*,p*) = 0. Now we use the same approach as above and linearize around 
the fixed point: 

q = p-p*,r = m-m*,q T =p T -p* (24) 
which leads to the following dynamical equations for the increments: 

q = f(p*+q,m*+r) = f(p*,m*) + ^-q+^-r = 

= + a-q + fi-r (25) 



r = d(P* + q,m* +r,p* +q T ) = 

i * * * A . 9g dg dg 
= 9(P ,m ,p ) + — ■ q+ — ■ r + — ■ q T = 
op dm op T 

= + "/-q + 5-r + e-q T (26) 
We can now assume solutions on the form: 

q(t) = qi e Xlt + q 2 e X2t 1 r{t) = r x e Xlt + r 2 e A2 * (27) 
and find after some lengthy algebra the eigenvalues: 

a + 5±((a-8f + 4/3(ee" AlT +7))^ (28) 



a + S±((a- S) 2 +A(3{ee- X2T +7))* (29) 



We note that these are transcendal equations in the eigenvalues Ai, A2 but that 
the time delay t specifically appears into the relations (for a further discussion, 
see previous Appendix). 
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In many biological systems with genetic feed back regulations, there are 
often more than two dynamical variables. We shall in the following sections 
discuss oscillatory states of the transcription factor NF-kB. As we show, that 
system can favourably be reduced to a three dimensional dynamical system in 
the variables NF-kB, its inhibitor IreB and the associated mRNA, IkBul Let 
us for simplicity write these variables as N, I, I m thus formally obtaining the 
following three dimensional dynamical system: 

N \ ( f(N,I m ,I) \ 
L = g(N,i m ,i) (30) 
/ J V HN,I m ,I) J 

The procedure to study this type of quations is exactly the same as for the two- 
dimensional system. From the stationary point N*,I m ,I* we linearize around 
it using the Jacobian matrix. 



T — ®<L dg dg_ /oi \ 

J ~ aN r5f„ ar \ 01 > 




In this case one obtains three eigenvalues Ai,A2,A3 which are either all three 
real, or one is real and the two others complex conjugates. Again, in the last 
case one may encounter a Hopf bifurcation where the fixed point N* , I* n , I* 
bifurcates into a limit cycle in the three-dimensional N,I m ,I space, causing 
oscillations in these variables. Notice, that there topologically is a big difference 
between oscillatory limit cycles in two and three dimensions. In two dimensions 
a limit cycle can 'enclose' trajectories, which is not the case in three dimensions 
where tracjectories may 'escape around' the limit cycles. An important theorem 
for this type of behavior is called the Poincare-Bendixson theorem, which states 
that if a trajectory is confined to a closed, bounded region and there are no 
fixed points in that region, then the trajectory eventually approaches a closed 
orbit (see e.g. [28]). 



C Stability analysis of delay systems 

C.l Linear delay systems 

Consider the simplest delayed system of the form 

dx 

— = ax(t) + bx(t-r). (32) 

A complete description of the solution is given in an appendix of ref. [65] . 
while we will just investigate the conditions where the system display sustained 
oscillations. The general solution of Eq. (|32[) is x = exp(Xt) which, inserted in t 
investigate the conditions where the system display sustained oscillations. The 
general solution of Eq. (|32|) gives 

A = a + 6e" AT . (33) 
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Defining A = /i + iui, the oscillatory case takes place when fi = 0. The solution 
is then 

iu> = a + b cos lot — ib sin lot (34) 

which gives 

(6 2 -a 2 ) 1/2 (35) 
arccos(—a/b) 

(6 2 _ a 2)l/2 • ( 36 ) 

A more general treatment (55] gives the conditions under which the system 
converges to a stationary solution, that is \a\ > \b\ for any r or if \a\ < \b\ for 
r < arccos(—a/b)/((b 2 — a 2 ) 1 / 2 ). Roughly speaking, the system can oscillate 
if the dominant part of the kernel is the delayed one and if the delay is large 
enough. 



LO = 
T = 



C.2 Absence of closed orbits 

If we label the vector forming the left-hand of Eq. JTJ) as x and that forming 
the side right-hand side as F, its divergence of F is 

V.F=g-*. + |-V (37) 

Since we exclude that each of the two molecules can activate itself, the partial 
derivatives df/dx and dg/dy are non-positive, and consequently V • F < 0. By 
virtue of Green's theorem, if the trajectory C were closed, 

> J V • idA = <j> x • n dl, (38) 

where n is the normal to the trajectory, and consequently x • n = everywhere. 
Since the circulation of a null vector is zero, this leads to a contradiction, and 
the trajectory cannot be closed. 



C.3 Stability analysis of the p53-mdm2 system 

Neaumtu and coworkers develop in ref. [JS] the complete stability analysis of 
the p53-mdm2 system with delay. First they show that for any choice of the 
parameters there is a unique stationary point for the rate equations. Consider 
the system obtained by linearization of Eq. ([5]) aroud such point 

(apio - b)p(t) - ap Q1 m(t) 

-dm(t) + cj w p(t - t) + c7oip(£ - r), (39) 

where pio and poi are the derivatives of pm(t) with respect to p and m in 
the stationary point, while 710 and 710 are the derivatives of the function p — 



dp 

di 
dm 

~dt 
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pm/ (kg + p — pm). Define p\ = b + d + apio, Po = db + adpia, qi = C701 and 
Qo = c7oi(6 + apio) — acpoi7io- Let A be the eigenvalues of the linearized system. 
It is possible to prove that if r > and p\ — q\p\ — 2p n q 2 > then there is a 
delay tq such that 



and consequently a Hopf bifurcation occurs at the stationary point. 

In the limit of large dissociation constant k, p 10 = P01 = 0. Consequently 
the condition for a Hopf bifurcation begins b 2 > 0, which is always satisfied 
(except for the trivial case b = 0). 

The critical delay tq is given by 




(40) 



To = — arcsin 
w V 



Pi^o 



+ arcsin 



qi^o 



(41) 



where wo is the solution of the equation 



4 . / 2 o 1 2\ 2 . 2 2 n 

u> + {-p 1 - 2p a + qi )uj +p -q =Q. 



(42) 
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